######################################################################
### Fig 2: Inconsistencies Between Heckscher-Ohlin Model and Actual Trade Flows
######################################################################
rm(list=ls())

library(foreign)
library(countrycode)
library(vcd)
library(png)

## setting the working directory to the current folder
setwd(getwd())

ternary <- function (x, scale = 1, dimnames = NULL,
                     dimnames_position = c("corner", "edge", "none"),
                     dimnames_color = "black", id = NULL, id_color = "black", 
                     id_just = c("center", "center"), coordinates = FALSE, grid = TRUE, 
                     grid_color = "gray", labels = c("inside", "outside", "none"), 
                     labels_color = "darkgray", border = "black", bg = "white", 
                     pch = 19, cex = 1, prop_size = FALSE, col = "red", main = "ternary plot", 
                     newpage = TRUE, pop = TRUE, ...) 
{
  labels <- match.arg(labels)
  if (grid == TRUE) 
    grid <- "dotted"
  if (coordinates) 
      id <- paste("(", round(x[, 1] * scale, 1), ",",
                  round(x[, 2] * scale, 1), ",", round(x[, 3] * scale, 1), ")", 
                sep = "")
  dimnames_position <- match.arg(dimnames_position)
  if (is.null(dimnames) && dimnames_position != "none") 
    dimnames <- colnames(x)
  if (is.logical(prop_size) && prop_size) 
    prop_size <- 3
  if (ncol(x) != 3) 
    stop("Need a matrix with 3 columns")
  if (any(x < 0)) 
    stop("X must be non-negative")
  s <- rowSums(x)
  if (any(s <= 0)) 
    stop("each row of X must have a positive sum")
  x <- x/s
  top <- sqrt(3)/2
  if (newpage) 
    grid.newpage()
  xlim <- c(-0.03, 1.03)
  ylim <- c(-1, top)
  pushViewport(viewport(width = unit(1, "snpc")))
  if (!is.null(main)) 
    grid.text(main, y = 0.9, gp = gpar(fontsize = 20, fontstyle = 1))
  pushViewport(viewport(width = 0.8, height = 0.8, xscale = xlim, 
                        yscale = ylim, name = "plot"))
  eps <- 0.01
  grid.polygon(c(0, 0.5, 1), c(0, top, 0), gp = gpar(fill = bg, 
                                             col = border), ...)
  if (dimnames_position == "corner") {
    grid.text(x = c(0, 1, 0.5), y = c(-0.04, -0.04, top + 
                                  0.07), label = dimnames, gp = gpar(fontsize = 13))
  }
  if (dimnames_position == "edge") {
    shift <- eps * if (labels == "outside") 
      8
    else 0
    grid.text(x = 0.25 - 2 * eps - shift, y = 0.5 * top + 
              shift, label = dimnames[2], rot = 60, gp = gpar(col = dimnames_color))
    grid.text(x = 0.75 + 3 * eps + shift, y = 0.5 * top + 
              shift, label = dimnames[1], rot = -60, gp = gpar(col = dimnames_color))
    grid.text(x = 0.5, y = -0.02 - shift, label = dimnames[3], 
              gp = gpar(col = dimnames_color))
  }
  if (is.character(grid)) 
    for (i in 1:4 * 0.2) {
      grid.lines(c(1 - i, (1 - i)/2), c(0, 1 - i) * top, 
                 gp = gpar(lty = grid, col = grid_color))
      grid.lines(c(1 - i, 1 - i + i/2), c(0, i) * top, 
                 gp = gpar(lty = grid, col = grid_color))
      grid.lines(c(i/2, 1 - i + i/2), c(i, i) * top, gp = gpar(lty = grid, 
                                                       col = grid_color))
      if (labels == "inside") {
        grid.text(x = (1 - i) * 3/4 - eps, y = (1 - i)/2 * 
                  top, label = i * scale, gp = gpar(col = labels_color), 
                  rot = 120)
        grid.text(x = 1 - i + i/4 + eps, y = i/2 * top - 
                  eps, label = (1 - i) * scale, gp = gpar(col = labels_color), 
                  rot = -120)
        grid.text(x = 0.5, y = i * top + eps, label = i * 
                  scale, gp = gpar(col = labels_color))
      }
      if (labels == "outside") {
        grid.text(x = (1 - i)/2 - 6 * eps, y = (1 - i) * 
                  top, label = (1 - i) * scale, gp = gpar(col = labels_color))
        grid.text(x = 1 - (1 - i)/2 + 3 * eps, y = (1 - 
                    i) * top + 5 * eps, label = i * scale, rot = -120, 
                  gp = gpar(col = labels_color))
        grid.text(x = i + eps, y = -0.05, label = (1 - 
                    i) * scale, vjust = 1, rot = 120, gp = gpar(col = labels_color))
      }
    }
  xp <- x[, 2] + x[, 3]/2
  yp <- x[, 3] * top
  size = unit(if (prop_size) 
    prop_size * (s/max(s))
  else cex, "lines")
  grid.points(xp, yp, pch = pch, gp = gpar(col = col), default.units = "snpc", 
              size = size, ...)
  if (!is.null(id)) 
    grid.text(x = xp, y = unit(yp - 0.015, "snpc") - 0.5 * 
              size, label = as.character(id), just = id_just, gp = gpar(col = id_color, 
                                                                cex = cex))
  if (pop) 
    popViewport(2)
  else upViewport(2)
}


### Imports 

## Testing the function

### 1991
result.imp <- read.csv("./1991_impHML.csv")
## remove outliers
lo.imp <- as.vector(quantile(result.imp$total, probs=c(.001,.999,1))[1])
med.imp <- as.vector(quantile(result.imp$total, probs=c(.001,.999,1))[2])
result.imp <- subset(result.imp, subset=(result.imp$total>=lo.imp & result.imp$total<=med.imp))

colnames(result.imp)[3] <- "Medium Wage States"
colnames(result.imp)[4] <- "Low Wage States"
colnames(result.imp)[5] <- "High Wage States"


### 2000
result.imp2000 <- read.csv("./2000_impHML.csv")
## remove outliers
lo.imp <- as.vector(quantile(result.imp2000$total, probs=c(.001,.999,1))[1])
med.imp <- as.vector(quantile(result.imp2000$total, probs=c(.001,.999,1))[2])
result.imp2000 <- subset(result.imp2000, subset=(result.imp2000$total>=lo.imp & result.imp2000$total<=med.imp))

colnames(result.imp2000)[3] <- "Medium Wage States"
colnames(result.imp2000)[4] <- "Low Wage States"
colnames(result.imp2000)[5] <- "High Wage States"

### 2009
result.imp2009 <- read.csv("./2009_impHML.csv")
## remove outliers
lo.imp <- as.vector(quantile(result.imp2009$total, probs=c(.001,.999,1))[1])
med.imp <- as.vector(quantile(result.imp2009$total, probs=c(.001,.999,1))[2])
result.imp2009 <- subset(result.imp2009, subset=(result.imp2009$total>=lo.imp & result.imp2009$total<=med.imp))

colnames(result.imp2009)[3] <- "Medium Wage States"
colnames(result.imp2009)[4] <- "Low Wage States"
colnames(result.imp2009)[5] <- "High Wage States"


png(file="./ternary_imports_final")
ternary(result.imp[,3:5], prop_size = TRUE, scale=TRUE,
        main="1991")
dev.off()

png(file="./ternary_imports_final2")
ternary(result.imp2000[,3:5], prop_size = TRUE, scale=TRUE,
        main="2000")
dev.off()

png(file="./ternary_imports_final3")
ternary(result.imp2009[,3:5], prop_size = TRUE, scale=TRUE,
        main="2009")
dev.off()


pdf(file="./Figure2a.pdf",
    width=18, height=6)

par(mar=c(0,4,0,0))
par(mfrow=c(1,3))
ima <- readPNG("./ternary_imports_final")
# Set up the plot area
par(mar=c(4,6,4,4))
par(omar=c(0,0,0,0))
plot(1:2, type='n', main="", xlab="", ylab="",
     xaxt='n', yaxt='n', bty='n')
# Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(ima, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
mtext(side=2, "Sources of US Imports", line=4, cex=2)
rect(1,1.04, 1.97,1.1,col = rgb(0.5,0.5,0.5,1/3),
     border=NA)

text(1.77, 1.44, labels="Heckscher-Ohlin", cex=2.1, col="gray20")
text(1.77, 1.39, labels="Prediction", cex=2.1, col="gray20")

arrows(1.8, 1.35, 1.5, 1.1, length = 0.25, angle = 30,
       code = 2, cex=.9, col="gray45")


ima2 <- readPNG("./ternary_imports_final2")
# Set up the plot area
par(omar=c(0,0,0,0))
plot(1:2, type='n', main="", xlab="", ylab="",
          xaxt='n', yaxt='n', bty='n')
#Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(ima2, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
rect(1,1.04, 1.97,1.1,col = rgb(0.5,0.5,0.5,1/3),
     border=NA)

ima3 <- readPNG("./ternary_imports_final3")
# Set up the plot area
par(omar=c(0,0,0,0))
plot(1:2, type='n', main="", xlab="", ylab="",
     xaxt='n', yaxt='n', bty='n')
# Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(ima3, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

rect(1,1.04, 1.97,1.1,col = rgb(0.5,0.5,0.5,1/3),
     border=NA)


dev.off()


### Exports 

## 1991
result.exp <- read.csv("./1991_expHML.csv")
## remove outliers
lo.exp <- as.vector(quantile(result.exp$total, probs=c(.001,.999,1))[1])
med.exp <- as.vector(quantile(result.exp$total, probs=c(.001,.999,1))[2])
result.exp <- subset(result.exp, subset=(result.exp$total>=lo.exp & result.exp$total<=med.exp))

colnames(result.exp)[3] <- "Medium Wage States"
colnames(result.exp)[4] <- "Low Wage States"
colnames(result.exp)[5] <- "High Wage States"


### 2000
result.exp2000 <- read.csv("./2000_expHML.csv")
## remove outliers
lo.exp <- as.vector(quantile(result.exp2000$total, probs=c(.001,.999,1))[1])
med.exp <- as.vector(quantile(result.exp2000$total, probs=c(.001,.999,1))[2])
result.exp2000 <- subset(result.exp2000, subset=(result.exp2000$total>=lo.exp & result.exp2000$total<=med.exp))

colnames(result.exp2000)[3] <- "Medium Wage States"
colnames(result.exp2000)[4] <- "Low Wage States"
colnames(result.exp2000)[5] <- "High Wage States"


### 2009
result.exp2009 <- read.csv("./2009_expHML.csv")
## remove outliers
lo.exp <- as.vector(quantile(result.exp2009$total, probs=c(.001,.999,1))[1])
med.exp <- as.vector(quantile(result.exp2009$total, probs=c(.001,.999,1))[2])
result.exp2009 <- subset(result.exp2009, subset=(result.exp2009$total>=lo.exp & result.exp2009$total<=med.exp))

colnames(result.exp2009)[3] <- "Medium Wage States"
colnames(result.exp2009)[4] <- "Low Wage States"
colnames(result.exp2009)[5] <- "High Wage States"

png(file="./ternary_exports_final")
ternary(result.exp[,3:5], prop_size = TRUE, scale=TRUE,
        main="1991")
dev.off()

png(file="./ternary_exports_final2")
ternary(result.exp2000[,3:5], prop_size = TRUE, scale=TRUE,
        main="2000")
dev.off()

png(file="./ternary_exports_final3")
ternary(result.exp2009[,3:5], prop_size = TRUE, scale=TRUE,
        main="2009")
dev.off()


pdf(file="./Figure2b.pdf",
    width=18, height=6)

par(mfrow=c(1,3))
par(mar=c(0,4,0,0))
par(omar=c(0,0,0,0))

ima <- readPNG("./ternary_exports_final")
# Set up the plot area
par(mar=c(4,6,4,4))
par(omar=c(0,0,0,0))
plot(1:2, type='n', main="", xlab="", ylab="",
     xaxt='n', yaxt='n', bty='n')
# Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(ima, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
mtext(side=2, "Destinations of US Exports", line=4, cex=2)
rect(1,1.04, 1.97,1.1,col = rgb(0.5,0.5,0.5,1/3),
     border=NA)

text(1.77, 1.44, labels="Heckscher-Ohlin", cex=2.1, col="gray20")
text(1.77, 1.39, labels="Prediction", cex=2.1, col="gray20")
arrows(1.8, 1.35, 1.5, 1.1, length = 0.25, angle = 30,
       code = 2, cex=.9, col="gray45")

ima2 <- readPNG("./ternary_exports_final2")
# Set up the plot area
par(omar=c(0,0,0,0))
plot(1:2, type='n', main="", xlab="", ylab="",
          xaxt='n', yaxt='n', bty='n')
# Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(ima2, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
rect(1,1.04, 1.97,1.1,col = rgb(0.5,0.5,0.5,1/3),
     border=NA)

ima3 <- readPNG("./ternary_exports_final3")
# Set up the plot area
par(omar=c(0,0,0,0))
plot(1:2, type='n', main="", xlab="", ylab="",
     xaxt='n', yaxt='n', bty='n')
# Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(ima3, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

rect(1,1.04, 1.97,1.1,col = rgb(0.5,0.5,0.5,1/3),
     border=NA)

dev.off()


